
conf = omni_load_conf()

myps,'./grs_paper/detection_flux_plot.eps',xsize=6.5,ysize=2.3,thick=1.5
!x.omargin = [-5,-6.5]
!y.omargin = [-0.5,-1.5]+0.25
multiplot,[2,1]

csf = 0.55

FOR jj=0,1 DO BEGIN
   
   s    = omni_read_cat(conf.cat,ncat)
   
   restore,'./local/'+conf.survey+'_velocities.sav',/ver
   
   ind = jj EQ 0 ? where(s.flux40 GT 0 AND s.glon GE 7.5 AND $
                         s.glon LE 194.,np) : $
         where(s.flux40 GT 0 AND s.glon GE 18 AND s.glon LE 56, np)
   
   v = v[ind]
   s = s[ind]
   
   
   logs = alog10(s.flux40)
   yarr = histogram(logs,loc=xarr,min=-2,bin=0.2,max=2,REVERSE=rind)
   
   nx = n_elements(xarr)
   dx = arrstep(xarr)
   print,dx
   detfg = dblarr(nx)
   detfd = dblarr(nx)
   detfa = dblarr(nx)
   
   FOR ii=0,nx-1 DO BEGIN       ; Loop through flux-density bins
      
      IF rind[ii] EQ rind[ii+1] THEN CONTINUE ; If nothing here, skip!
      
      inds = rind[rind[ii]:rind[ii+1]-1] ; Extract the inds of the relevant objs
      nr = n_elements(inds)              ; Count
      
      vel = v[inds]             ; Make easy to type
      
      gi = where(vel.grs.flag EQ 1 OR vel.grs.flag EQ 2 OR $ ; Indices for GRS
                 vel.grs.flag EQ 5 OR vel.grs.flag EQ 6, ng) ; objects
      
      di  = where(vel.vlsr GE -500, nd)  ; Indices for dense-gas objects
      din = where(vel.vlsr LT -500, ndn) ; Indices for BAD dense-gas objects
      
      CASE 1 OF                 ; Concatenate these for the "ALL INDICES"
         ng NE 0 AND nd NE 0: ai = [gi,di]
         ng NE 0: ai = gi
         nd NE 0: ai = di
         ELSE: ai = !null
      ENDCASE
      
      IF n_elements(ai) GT 0 THEN BEGIN ; If there are any...
         ai = ai[uniq(ai,sort(ai))]
         ni = missing(ai,lindgen(nr),nn)   
      ENDIF ELSE BEGIN          ; Else, there are none
         ni = lindgen(nr)
         nn = nr
      ENDELSE
      
      na = n_elements(ai)
      
      ni = missing(ai,lindgen(nr),nn)   
      
      detfg[ii] = double(ng) / nr
      detfd[ii] = double(nd) / nr
      detfa[ii] = double(na) / nr
      
      IF ng NE 0 THEN glon = m4_stat(vel.l) ELSE glon = [0,0,0,0]
      IF ng NE 0 THEN glat = m4_stat(vel.b) ELSE glat = [0,0,0,0]
      
      print,'------------'
      print,xarr[ii]+dx/2.,detfa[ii],nr
      print,glon
      print,glat
      ;; IF xarr[ii]+dx/2. GE 0.6 THEN BEGIN
      ;;    print,'************'
      ;;    print,nd,ndn,nr
      ;;    FOR kk=0,ndn-1 DO BEGIN
      ;;       print,vel[din[kk]].l,vel[din[kk]].b,$
      ;;             vel[din[kk]].cnum,vel[din[kk]].mol[0].nspec
      ;;       quick_view_bgps_hht,vel[din[kk]].l,vel[din[kk]].b
      ;;       res = get_kbrd()
      ;;    ENDFOR
      ;; ENDIF
      
   ENDFOR                       ; End of loop over flux density bins
   
   
   cgPlot,xarr+dx/2.,detfd,color='black',linestyle=3,$
          charsize=1.0*csf,psym=10,yr=[0,1.05],xr=[-1.8,1.6],thick=2.5,$
          ytit=jj EQ 0 ? 'V!dLSR!n Assignment Fraction / CDF' : '',xst=5,/nodata
 
   ;; overplot the flux-density CDF
   xp = alog10(s.flux40)
   yp = ccdf(xp,/cdf)
   cgOplot,xp,yp,color='blu4',linestyle=2
   
   cgOplot,xarr+dx/2.,detfd,color='blk5',linestyle=0,psym=10,thick=2.5
   IF jj EQ 1 THEN BEGIN
      cgOplot,xarr+dx/2.,detfg,color='black',psym=10,thick=1.5,linestyle=3
      cgOplot,xarr+dx/2.,detfa,psym=10,thick=1.5
      forprint,xarr+dx/2.,detfa
   ENDIF
   
   cgAxis,xaxis=0,charsize=1.0*csf,xtit='Flux Density  [Jy]',/xlog,/xst,/save
   cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
   cgAxis,yaxis=0,/yst,ytickformat='blank_axis'
   cgAxis,yaxis=1,/yst,ytickformat='blank_axis'
   
   leq = cgSymbol('leq')
   deg = cgSymbol('deg')
   
   CASE jj OF
      0 : BEGIN
         al_legend,/top,/left,box=0,charsize=1.0*csf,$
                   ['7.5'+deg+' '+leq+' l '+leq+' 194'+deg]
      END
      1: BEGIN
         al_legend,/bottom,/right,box=0,charsize=0.9*csf,linsize=2.0,$
                   linestyle=[0,3,0,2],$
                   thick=[5,3,3,3]/2.,color=['blk5','black','black','blu4'],$
                   ['Dense Gas Tracers','GRS !u13!nCO',$
                    'Combined','Flux Density CDF'],/clear,spacing=0.8
         al_legend,/top,/left,box=0,charsize=1.0*csf,$
                   ['18'+deg+' '+leq+' l '+leq+' 56'+deg]
      END
   ENDCASE 
      
   multiplot
   
ENDFOR

myps,/done,/mp

END
